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ABSTRACT 


Model reduction is an important practical problem in the control of flexible spacecraft, 
and a considerable amount of work has been carried out on this topic. Two of the best- 
known methods developed are modal truncation and internal balancing. Modal 
truncation is simple to implement but can give poor results when the structure possesses 
clustered natural frequencies, as often occurs in practice. Balancing avoids this problem 
but has the disadvantages of high computational cost, possible numerical sensitivity 
problems, and no physical interpretation for the resulting balanced "modes". 

The purpose of this work is to examine the performance of the subsystem balancing 
technique developed by the investigator when tested on a realistic flexible space 
structure, in this case a model of the Permanently Manned Configuration (PMC) of Space 
Station Freedom. This method retains the desirable properties of standard balancing 
while overcoming the three difficulties listed above. It achieves this by first 
decomposing the structural model into subsystems of highly correlated modes. Each 
subsystem is approximately uncorrelated from all others, so balancing them separately 
and then combining yields comparable results to balancing the entire structure directly. 
The operation count reduction obtained by the new technique is considerable: a factor of 
roughly if the system decomposes into r equal subsystems. Numerical accuracy is also 
improved significantly, as the matrices being operated on are of reduced dimension, and 
the modes of the reduced-order model now have a clear physical interpretation; they are, 
to first order, linear combinations of repeated-frequency modes. 
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INTRODUCTION 


Model reduction is a very important practical problem related to the control of flexible 
space structures (FSS), and a considerable amount of work has been earned out on this 
topic. Well-known methods include modal truncation [1], based either on the natural 
frequencies of the structure or its modal costs, and balancing [2] of the entire structure 
and then truncation to retain a dominant model for it An advantage of the balancing 
approach is that it typically yields a more accurate reduced-order model than does simple 
modal truncation. This is particularly true when the structure possesses clustered natural 
frequencies, as is often the case for realistic flexible space structures. However, the 
disadvantages of balancing are its high computational cost, possible numerical sensitivity 
problems resulting from the large matrices being operated on, and the difficulty involved 
in providing a physical interpretation for the resulting balanced "modes". 

The purpose of this paper is to investigate the practical performance of the alternative 
subsystem balancing technique when tested on a realistic flexible space structure. This 
method, introduced in [3] [4] and further developed in [5], retains the desirable properties 
of standard balancing while overcoming the three difficulties listed above. This is 
achieved by first decomposing the structural model into subsystems of highly correlated 
modes, based on the modal correlation coefficients derived in [4] from the controllability 
and observability Grammian matrices [6] of the structure. Each subsystem is 
approximately uncorrelated from all others, so balancing each separately and 
concatenating the dominant reduced-order models obtained yields roughly the same result 
as balancing the entire structure directly. The computational cost reduction produced by 
this block-by-block technique is considerable: an operation count reduction by a factor of 
roughly y.t if the system decomposes into r equal subsystems. The numerical accuracy 
of die resulting reduced-order model is also improved considerably, as the matrices being 
operated on are of reduced dimension; this avoids the numerical conditioning problems 
noted in [8][9] for standard balancing. Furthermore, the modes of the reduced model do 
now permit a clear physical interpretation. This is a consequence of the fact that each 
correlated subsystem must necessarily only include modes with close natural frequencies. 
The balanced modes of each subsystem are therefore, to first order, linear combinations of 
repeated-frequency modes, and so can themselves be taken as an equally valid set of 
physical modes. Balancing the entire structure, on the other hand, combines modes of 
widely differing frequencies, making interpretation difficult 

The results obtained using the software described in this report are for the Permanently 
Manned Configuration (PMC) of Space Station Freedom. Two different "stick models" 
[11] for this vehicle were studied, for two choices of solar array and radiator orientations. 
In both cases, the initial 202-mode flexible body models could be reduced to models with 
between 20 and 30 modes with very little loss of accuracy. 

THEORETICAL BACKGROUND 

Consider an n-mode model for the structural dynamics of a modally damped, non- 
gyroscopic, non-circulatory FSS with m actuators and p sensors, not necessarily 
collocated. This model can be written in modal form [1] as 
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( 1 ) 


ii + diagil^Q),)^ + diag(Q)f) n = flu, 
y = C r f|+C rf ti, 


where t| is the vector of modal coordinates, u that of applied actuator inputs and y that of 
sensor outputs, and a>, and £, are the natural frequency and damping ratio of the i*h 
mode, respectively. For the typical FSS [7], the {Q are quite low (e.g. 0.5 %), and the 

{o>,} occur in clusters of repeated, or nearly repeated, frequencies as a result of structural 
symmetry. 

Defining the state vector x = (Th,a) l rj„“-,rj a ,a a t) l ,) T for this structure yields the state 
space representation x = Ax + flu,y = Cx, where A = blkdiag(A) it B = (B[,—,I%) T and 
C = (CJ,— ,CJ, with 


4 = 



?)«-(*) 


and Cj = (c rt c s l(o i )-. 


( 2 ) 


b, is the i* row of fl, and c ri and c s are the i* columns ofC, and Q, respectively. 
The problem studied here is that of obtaining a reduced-order model 


y r = c,x r 

for this structure for which the normalized output error 


(3) 


= Jly(0 - yr(')£<*/J|y(o£<* 


(4) 


is acceptably sm al l . Of course, the size of S will depend on the order, n r , chosen for the 
reduced model. A good model reduction procedure should ideally provide information 
allowing an intelligent choice for nr to be made so as to achieve a specified S value. 

Two techniques for model reduction that have been extensively studied are those of 
modal truncation and internal balancing. The new method im plemente d in this report, 
subsystem balancing, can be regarded as an intermediate case between the two established 
techniques. Model reduction by subsystem balancing proceeds by first dividing the given 
structure into subsystems of highly correlated modes. Each subsystem is then balanced 
independently, and a reduced-order model for it generated by deleting all balanced s ta t e s 
corresponding to Hankel singular values [2] below some specified threshold. Note that 
the singular value weighting described in [10] could be applied, if desired, without 
changing the argument in any way. Similarly, frequency w eightin g of the Hankel 
singular values can easily be incorporated to deal with input signals which have a known 
frequency spectrum. This is actually done in the present application, where the inputs are 
steps (representing thruster firings) rather than the impulses c lassic ally considered in 


22-4 



model reduction problems. The resulting reduced-order subsystem models so obtained 
are then combined to yield a dominant, approximately balanced, reduced-order model for 
the full system. 


USER INTERFACE 

This section describes the user interface to the model reduction package which was 
developed as part of this contract This software consists of a library of Matlab m- 
functions, with mrmain calling all the other functions internally. The package is installed 
on the Sun SparcStation 2 deimos in the Integrated Analysis Laboratory in Building 16, 
and hat also been produced in a Macintosh version. The documentation that follows 
details the user interface for mrmain; listings of this function, together with the second- 
level functions it calls, are given as an appendix. All functions have extensive in-line 
documentation, facilitating future use and/or modification. 

Input arguments 

ortv The natural frequencies (rad/s) of the structure, input as either a row or column 
vector. Any rigid-body modes must precede the flexible modes and be represented by 
hard zero frequencies. 

phia : The influence matrix, in mass-normalized coordinates, corresponding to the 
specified actuator locations. If the structure has n modes and m actuators, phia will be an 
(n x m) matrix. 

phis: Similar to phia, but for sensor stations or positions of outputs of interest (e.g. solar 
array tips). 

User responses 

Output the time taken for each step?: The time required for each matrix decomposition, 
etc., is output to the screen if requested. This allows the progress of the model reduction 
procedure to be monitored, as well as giving an indication of which steps are the most 
computationally intensive. 

Vectorize? (Faster, but requires more storage): In Matlab, for loops are typically an 
order of magnitude slower to excute than the equivalent "vectorized" operation. For 
instance, s*0; for t=1:n, s=s+x(i); end; runs considerably slower than does s=x*ones(n,1). If 
vectorization is requested, computation of the system Grammian matrices and correlation 
coefficients is put into the form of vector-matrix operations rather than loops; this is 
indeed considerably faster, but requires some additional temporary storage arrays. 

Structural damping ratio, % (default is 0.5%): The specified damping ratio is applied 
uniformly to all flexible modes of the full structural model. 

Print frequencies in Hz?: The mean frequency of each subsystem can be output in either 
rad/s or Hz, as desired. 

Desired controllability threshold?: This threshold value is used to determine which 
modes are correlated in a controllability sense. The system is then broken down into 
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disjoint sets of modes (subsystems), where modes with a controllability correlation 
coefficient greater than the specified threshold are deemed to be correlated. Taking a 
threshold value of 0 implies that all modes are considered correlated, i.e. the method 
reduces to standard balancing. Conversely, a threshold value of 1 implies that no modes 
are taken together this is modal truncation. Intermediate values allow the dimensions of 
the resulting subsystems to be specified to a large extent; reducing the threshold reduces 
the number of subsystems, so increasing their dimension. 

Desired overall threshold ?: This threshold is used in a similar fashion to the 
controllability threshold, but both controllability and observability are now taken into 
account This yields the final subsystem distribution output by the program (in modemap) 
and used to obtain the reduced-order model. 

Compare step responses ?: If requested, the step responses of the full and reduced-order 
models are computed, plotted, and the relative differences (Le. reduced-order model error) 
output for each input-output channel. 

Desired truncation measure?: Two types of measure can be used to define the number of 
modes retained in the reduced-order model. If a positive integer is input, this is taken to 
be precisely the desired reduced-order model. On the other hand, if a real number in the 
interval [0, 1) is input, this is taken to be the desired relative error in the reduced-order 
model step response, and the model order required to achieve this is computed. (Note 
that this later option is only an approximation, and should only be used as such.) 

Output arguments 

am, bm, cm : The reduced-order state-space model obtained. 

modemap: This matrix specifies which physical modes are grouped into which 
subsystems in the decomposition based on overall correlation coefficients. The i^ 1 
column of modemap lists the modes making up the i* of these subsystems. 

SUMMARY OF RESULTS 

Results will now be provided which illustrate the behavior of the subsystem balancing 
technique when applied to a structural model [11] of the Permanently M ame d 
Configuration (PMC) of Space Station Freedom. This structure pos se ss e s light damping 
(estimated to be 0.5% of critical), and a large number of closely-spaced vibration modes 
(202 flexible modes below 10 Hz). Two configurations of the PMC were investigated: in 
the first, the solar arrays are in the station yz-plane (a=/J = 0) and the radiators in the 
xy-plane (y = 0); in the second, the arrays are in the xy-plane and the radiators in the xz- 

plane (a = y = 90*, fi = 0). The inputs to these models are the 12 Reaction Control 
System (RCS) thrusters, Le. the port/starboard and upper/lower x, y and z jets. The 
measured outputs are the 3 angular rates sensed by the rate gyros on the station avionics 
pallet (The movements at other positions of interest for instanc e the solar array tips, 
could also be considered if desired; the method remains exactly the same.) 

A first point to examine is the efficiency of subsystem b alancin g as compared with that of 
standard balancing. Matlab function obalreal in the Robust Control Toolbox is a reliable 
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implementation of Moore’s balancing algorithm; applying this to the PMC models 
considered requires approximately 3 hours on a SparcStation 2. By contrast, the 
subsystem balancing implementation provided by mimain requires approximately 3 
minutes. Furthermore, the bulk of the operations in subsystem balancing are order(n 2 ), 
due in part to the use of closed-form Grammians [6], whereas standard balancing is 
order(n 3 ). The efficiency advantages of the new approach will therefore become only 
more pronounced as larger systems are examined. 

The role of the threshold coefficient in determining subsystem dimensions can be seen 
from the following table. The first column gives various choices for the controllability 
correlation threshold parameter, and columns 2 and 3 show the resulting maximum 
subsystem dimensions for the two PMC configurations studied (for input axis port upper 
x). It can be seen that these dimensions do indeed decrease as the threshold increases, as 
expected. Also, both systems exhibit broadly similar behavior. It can be noted that the 
evolution of subsystem dimensions is fairly discontinuous: for instance, large changes 
occur for thresholds between 0.10 and 0.15, whereas there are hardly any differences 
between 0.30 and 0.45. A consequence of this is that it is not always possible to find a 
threshold value which will yield a particular maximum subsystem order. However, it is 
possible to obtain a good working value which gives a totally acceptable subsystem 
partition. For the system studied here, a maximum subsystem dimension of about 30 
leads to about 36-38 individual subsystems (some of which consist of single modes), a 
good balance; threshold values giving this distribution were chosen as nominal. Using 
these thresholds, the original 202-mode flexible models were found to be reducible to 
models with only 20 to 30 modes without introducing significant errors into the resulting 
step responses. 


TABLE 1. - MAXIMUM SUBSYSTEM DIMENSIONS VERSUS THRESHOLD 


Threshold 

Max dim, a = 0' 

Max dim, a = 90* 

a=90\ C=l% 

0.00 

202 

202 

202 

0.05 

165 

165 

199 

0.10 

146 

131 

165 

0.15 

76 

87 

131 

0.20 

55 

51 

118 

0.25 

55 

31 

118 

0.30 

30 

20 

87 

0.35 

30 

20 

76 

0.40 

17 

20 

31 

0.45 

17 

18 

31 


The fourth column of the table illustrates the effect of damping on model reduction. 
Damping of flexible structures is a very difficult quantity to model, so there is 
considerable uncertainty in the damping levels to be chosen for Space Station Freedom. 

If a value of 1% of critical is used instead of the previous "nominal” lavel of 0.5%, the 
consequent broadening of the peaks of each mode increases the coupling between modes, 
so increasing the subsystem dimensions. This does not pose any problem, however; 
increasing the threshold value to 0.4 will again allow the desired dimensions to be 
obtained. 
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MATLAB PROGRAM LISTINGS 


function 

[am t bm,cm 4 nodcnu^>]=*mnnain(om,phU,phis); 

% 

% An M-function to perform model reduction based 
on 

% subsystem balancing for a uniformly-damped 
% flexible structure with rate outputs 
% 

% All other functions in this package are called 
% by this main routine. 

% 

% Arguments: 

% 

% In: Frequency vector om (rad/s), with any 
% rigid-body modes (om(i)=0) leading; 

% influence matrices phia (actuators) 

% and phis (sensors), one row per mode 
% 

% Out* Reduced-order state-space model {am,bm,cm} 
% of flexible-body dynamics; 

% the i-tfa column of modemap lists those modes 
% making up the i-th unreduced subsystem 
% 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n=max(size(om)); % Number of modes considered 
% 

% Strip off any rigid-body modes 
% 

nflex=sum(sign(om)); ruigid=n-nflex; 

om=om(nrigid+l :n); 

b=phia(nrigid+l:n,:); 

c=phis(nrigid+l:n,:)'; 

n=nflex; 

% 

timeout=0; 

stime=input( 1 Output the time taken for each step? \ 

V); 

if stime -* y timeout® 1; end; % Time o/p wanted 
% 

vect=0; 

svect=input( 'Vectorize? (Faster, but requires more 
storage) \ 's'); 

if svect *= y vect^l; end; % Vcctorizatioo wanted 
% 

% Enter specified damping ratio 
% 

ze=inputfStructural damping ratio, % (default is 
0.5%) *); 

if isempty(ze) = 1 ze= 03; end; ze=ze/100; 

% 

radhz=l; 

shz=inputCPrint frequencies in Hz? \ 's'); 

if shz ** *y' radhz=l/(2*pi); end; % Output format 

% 

% First compute controllability correlation coeffs 


% 

tO=clock; 

ro=cccalc(om,ze,b,vect); 
if timeout -= 0 tl=etime(clock,tO); 
s=r Ro-c completed after num2str(tl), * s']; 
dispC y,&*p(s)\ 
end; 

% 

% Next, find observability Grammian for entire sys 
% 

t0=clock; 

wo=cfgram(om,ze,c\- 1 ,vect); 
if timeout -= 0 tl=etime(clock,t0); 
s=[* Wo completed after ’, num2str(tl), * s']; 
dispC *); disp(s); 
end; 

% 

% Input desired controllability corrern threshold 
% 

dispC y 

dispC Threshold values should lie between 0 and 1;') 
dispC lower values give fewer (but larger) 
subsystems.*) 

dispC Enter a negative value when finished.') 

% 

xdum=l; % Dummy: gives inelegant indefinite loop 
while xdum > 0 

fn min acinputfDcs fared controllability threshold? '); 
if romin >■ 1 romin=l -eps,end; 
if romin >* 0 

% 

% Determine subsystems of correlated modes 
% 

tD=clock; 

[isort nsub]=subsys(ro jomin); 
if timeout ~= 0 tl=etime(clock,tO); 

*=[' Internal decomposition took ’, num2str(tl), ’ 

dispC *); disp(s); 
end; 

kmax=max(size(nsub)) ; % Num of subsystems 
«-[' Yields int2str(lanax), * subsystems; 
maximum size 1; 

e=[s, int2str(max(nsub)X minimum 
int2str(min(nsub))]; 
dispC '); disp(s); 

% 

else xdum = -l; 
end; 

% 

cad; 

% 

% Set up index to reorder Wo (agrees with isort) 

% 

iwsoit(2:2:2*n)=2*isort; 
iwsort(l :2:2*n-l )=2*isoit-ones( 1 ,n); 

% 
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% Operate on each subsystem in turn 

% 

iUl; 

nmaxzmax(nsub); % Greatest subsystem order 
wctoKl; 
tO-clock; 
for k*l±max 
nsubk*nsub(k); 

i2*il+nsubk-l; ivect*isort(il:i2); 
iwvecbriwsort(2*il-l:2*i2); 

% 

% First find its controllability Grammian 

% 


% 

% 

% 


wck*cfgram(om(ivect)je,b(ivect f :),l l vect); 

... then find its singular value decomposition 


[uk.sk, vk]asvd(wck); 
wck=uk*sqrt(ik); 

% 

% Finally, apply to correct row & col blocks of Wo 
% 


wo(iwvect,:)*rwck , ^wo(iwvect,:); 

wo(:,iwvect>*wo(:4Wvect)*wck; 


wctot*[wctot [wck; zcros(2+(nmax- 
nsubkX2*nsubk)]]; 
il*i2+l; 
end; 

if timeout -« 0 tl*etime(clock,tO); 
i*r Wbar completed after \ num2str(tl), * *1; 
d“PC *); 
end; 

% 

% Find matrix of overall correlation coefficients 

% 

tO=clock; 

ro=K>ccalc(wo,vect); 

if timeout 0 1 l*ctime(ckxdc,tO); 

**T Ro-o completed after num2str(tl) 1 * si; 

dispf 1; disp(s); 

end; 

% 

% Now try various overall threshold values 
% 

xdum «1; % Dummy variable again 
while xdum > 0 

romin2-inpqtCDesired overall threshold? '); 
if romin2 >* 1 romin2^ 1 -cps,end; 
if romin2 >« 0 

% 

% Determine new subsystems of correlated inodes 
% 

tO-clock; 

[isort2 nsub 2 }*subsys(rD 4 omin 2 ); 
if timeout 0 tWtime(clock,tO); 

**T Final decomposition took ', num2str(tl), 1 

* 1 ; 

dispf ’); disp(s); 
end; 

kmax2 »max(«ize(nsub2)); %Num of subsystems 


•“[’ Yields int2str(kmax2), * subsystems; 
maximum size 1; 

*=[*, int2str(max(nsub2)), minimum ’ t 
int2str(min(nsub2»]; 
dispC ’); disp(s); 

% 

else xdum = -1; 
end; 

% 

end; 

% 

% Now that the final subsystems are defined, 
% need to define corresponding Wo ordering 


iwsort2(2:2;2*n)=2*isort2; 
iwsort2(l :2:2*n-l )«2*isort2-ones(l uj); 

% 

% Find balancing transformation for each subsystem; 
% store it and the weighted Hankel singular values 

il-1; 

nmax=max(nsub2); % Greatest subsystem order 

omb«r=zeros( 1 Janax2); 

modemap=zeros(nmaxJanax2); 

bsv2tstQ; 

ttotsQ; 

tO*clock; 

% 

for k«l±tnax2 
nsubk*nsub2(k); 

i2*il+nsubk-l; ivect=isort2(il:i2); 

iwvect»iwsort2(2*il-l :2*i2); 

m od em a p ( 1 :nsubkjc>rivect , + nrigid*ones(nsubk, 1); 

% 

% Solve this sym eigen -/singular value problem 


% 

% 

% 


(tkJisv2k,vk]=svd(wo(iwvect > iwvect)); 

hsv2k*diag(hsv2k); 

[hsv2k ihsv]=sort(-hsv2k); hsv2k=-hsv2k; 
tk*tk(:ohsv); % Soil in descending order 
ttoHdot [tk; zeros(2*(nmax-nsubk),2*nsubk)]]; 


Perform ad hoc step response frequency weighting 


ombar(k)vsum(om(ivect)ynsu bk; 
hsv2ksfasv2k/(ombar(k)*ombar(k)); 
hsv2t»[hsv2t; hsv2k]; 

% 


1M2+1; 

«od; 

hsv2sum*sum(hsv2t); 

[hsv2s ihsv]*sort(-hsv2t); hsv2s*-hsv2s; 
if timeout 0 tl=etime(ck>ck,tO); 

**T Hankel singular values took \ num2str(tl), * si; 

dispC *); <fisp(s); 

end; 

% 

% Compute the balanced full-order system 

% 

clear wo; % Need the space for am (2nx2n here!) 
tO-clock; 
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[am bm cmj^modal(oin^e,b,c); 
if timeout ~= 0 tl«etime(clock,tO); 

*=[’ Full modal model took \ num2str(tl), ‘ >1; 

dispC y, di*p(s); 

end; 

% 

% Now apply Wc similarity transformation, by blocks 
% 

tO=clock; 

[am bm cm]=blkmult(am,bm,cm,wctot4Wsort4isub); 
if timeout 0 tl*etime(ck>ck,tO); 
s=C Wc similarity took \ num2str(tl), ’ si; 
dispC l;di sp(s); 
end; 

% 

% Then apply T similarity transformation, by blocks 
% 

tO=clock; 

[am bm cm]~blkmult(am,bm,cm v ttot4wsort24ksub2); 
if timeout ~= 0 tl«etime(ck>ck,tQ); 
s*r T similarity took \ num2str(tl), ’ si; 
dispC 1; disp(s); 
end; 

% 

% Set up for step response calculations, if wanted 
% 

dispc y 

sstep=input( 'Compare step responses? \ 's’); 
if sstep = y 

t=4/(ze*min(om)); % Tune for decay 
t=100*n>und(t/100); % Round to nearest hundred 
dt=t/100; t~[dt:dt:t]; % 100 points 
% 

m=size(b); m=m(2); % Number of inputs 
p=size(c); p=p(l); % Number of outputs 
% 

% Compute and store step responses of full system 

yMJ; 

tO=clock; 
for iu=l:m 

yf=[yf modstep(oin,ze,b # c,iu,t)] ; 
end; 

% 

if timeout -* 0 tl=etime(clock,tO); 
s=f Full step response took \ num2str(tl), ' si; 
dispC *); disp(s); 
end; 
end; 

% 

% Try various different truncation measures 
% 

dispC 1; 

dispC Enter either the desired number of modes 
(integer >0)') 

dispC or the acceptable approx relative output error (< 

DO 

dispC enter a negative quantity when finished.*) 

% 

xdum = 1; % Dummy variable again 
while xdum > 0 

cutoff=input('Desired truncation measure? 1; 


if isempty(cutoff) = 1 cutoff=l; end; % Safety 
t0=clock; 

if cutoff <0 xdum = -1; 
else 

if cutoff >= 1 nrom=min(cutoff,n); % # modes 
else % Find num modes from desired rel err 
abserr2=cutoff*cutoff*hsv2sum; 
test=0; i=n; 
while test <= abserr2 

test=test+hsv2s(2*i)+hsv2s(2*i- 1 ); 

W-l; 

end; 

nromsi+1; 

end; 

% 

% Find number of modes kept from each subsys 
% 

hsv2min=hsv2s(2*nrom); 
nsubr=zeros( 1 Janax2); 
ioff=0; 

for k=l:kmax2 
fori=l:nsub2(k) 

if abs(hsv2t(2*(ioff+i))) >= hsv2min 
nsubr(k)=nsubr(k)+ 1 ; 
end; 
end; 

ioff=ioff+nsub2(k); 

end; 

% 

% Finally, find truncation index to give ROM 

% 

il-1; 

iavect^G; 

% 

fork=ldcmax2 

i2=il+2*nsubr(k)-l ; 
iavect=[iavect iwsort2(il:i2)]; 
il«il+2*nsub2(k); 
end; 

% 

% Alter cm so as to give zero steady-state error 
% 

delcm=0*cm; 

g=am(iavect4avect)\bm(iavect,:); 
x=*-cm(:,iavect)*g/(g , *g); 
delcm(: 4 avect)*x*g’; 

% 

if timeout 0 tl*etime(ck>ck,tO); 
s»f Model dimension found after \ 
num2str(tlX 4 si; 

dispC l;disp(s); 
end; 

% 

% Output subsystem information 
% 

dispC 1; 

dispC Subsystem Number Mean freq’); 
s =£* dimension retained ’; 
if shz — y s=[s, ’ (Hz)l; else s=[s, *(rad/s)l; end 
disp(s); 

% 
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for k-l±max2 

if nsub2(k) >= 10 spad-' *; cod; 
if nsub2(k) >= 100 spad-"; end; 

•-T ’• q»d, int2str(nsub2(k))] ; 

% 

sp(d-' *; 

if nsubrflc) >■= 10 spad-’ cad; 
if nsubr(k) >- 100 (pad-"; cod; 

H*. ' \ (pad, int2str{nsubr(k))] ; 

% 

^ad- f 

prombar-radhz*ombar(k); 
if prombar >- 10 (pads’ end; 
if prombar >* 100 spad- 1 ’; cad; 
if prombar >- 1000 spad-"; cad; 

»*[»» # \ (pad, num2itr(prombar)]; 

dM(); 

ead; 

% 

% Compare step responses if requested 
% 

if (step — y 
for iu-l:m 
tO-clock; 

amp-am(iavect4avect); 
bmp=bm(iavect,:); 
cmp-cm(:4avect)^k»n(:4svect); 
yr*4>Dcstep(amp l bmp,anp,iu t t); 
if timeout — 0 tl-etime(clock,tO); 

®*T ROM step response took num2str(tl), 

dispC *); di*p((); 

rod; 

for k>-l:p 

iyf-(iu-l)*p+io; % Access correct 

response 

J^°*Cttyf(:4yf) yr(:4o)]); 

(*TStep responses, nrom - \ 
int2str(nrom)]; 

A». *; input *, int2str(iu)]; 
t-Cs, 1 , output \ int2str(k>)]; 

xlabclfTime (s) r ); 
ylabclfOutputt’); 


PMt,(yf(:4yf>-yr(:4o))); 

•*fStep error, mom « \ int2str(nrom)]; 
»-{*, *; input \ int2str(iu)l; 
c(s, f , output \ int2str(io)]; 
title(s); grid; 

xlabel(Time (a)’); 
yUbelCEnor'); 

pause; 

y e rr -o o rm (yf(:4yf)- 

yr(:,»)yaorm(yf(;4yf)); 

**f 2-norm relative output error of \ 
num2(tr(yerr)]; 

dispC *); disp(i); 


end; 

end; 

end; 

end; 

ead; 

% 

% Finally, store the chosen reduced-order model 
% 

am-am(iavect4avect); 

bm-bm(iavect,:); 

cm=cm(: 4avect>+delcm(: 4avect); 

* 

function ro-cccalc(om^e,b t vect); 

% An M-function to construct the controllability 
% correlation coefficients of a uniformly- 
% damped flexible structure 
% 

% The flag vect — 0 for Matlab-style vectorized 
% operations: faster, but requires extra temp store 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n— max(size(om)); % Number of modes considered 
% 

if vect — 0 
% 

% "Standard" loop over matrix locations 
% 

ro— eye(n); % Initialization (saves time) 

betad-zeros(n,l); % " " 

% 

% First compute contribution from B 
% 

for i-lm 
forj-i+lm 
ro(io>b(i,:)*b(j,:y; 

end; 

beUd(i)=*Mx(b(|,:)*b{v)’. eps); 

end 

% 

for i-lm 
forj-i+lm 

ro(ij)-abs(ro(ij)ysqrt{betad(i)*betad(j)); 

end; 

end 

% 

% Now add the frequency effects (Frob-aonn version) 
% 

for i-lm 
forj-i+lm 
g=om(jyom(i); 
temp-8*ze*xe*g; 
num-(g- l)*2*(g A 2+l); 
num*sqrt(temp*((g*teinp>+num)); 
d^Ks+WCg-l )^24<temp/2)); 


ro(ij)»=(num/den)*ro(ij); 
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ro(j,i)*ro(ij); 

end; 

end 

else 

% 

% Calculations in Matlab-style vectorized form: 

% first compute contribution from B 
% 

ro=b*b’; 

ro=abs(ro); 

% 

% Avoid singularities caused by zero entries in b f c 
% 

betad=max(diag(ro), eps*ones(n,l)); 
ro=ro-diag(diag(ro)); % Zero all diagonals 
ro=ro 4 diag(betad); % Put back diag, or eps 
% 

betads=sqrt(diag(betad)); 

ro=betads\abs(roVbetads; 

% 

% Now add the frequency effects (Frob-nonn version) 
% 

temp=8*ze*ze; 

% 

som=size(om); 

if som(l) >1 % om entered as a column vector 

g=(om .'X-l^om'; 

else % om entered as a row vector 

gsKom.^-ljy^om; 
end 
% 

num=((g-ooes(n)). A 2).*((g. A 2)+ooes(n)); 
num=sqrt((teinp*g).*((temp*(g. A 2))+num)); 
den^g+oi^n)).*(((g-ones(n)) A 2)4<temp*g/2)); 
rc*=(num7den).*ro; 
end 
# 


function w==cfgram(oDMe,b r cobs t vect); 

% 

% An M-function to compute the closed-form 
% Grammians of a unifonnly-damped flexible 
% structure with rate measurements 
% 

% The flag cobs “ l for controllability, 

% -1 for observability. 

% 

% The flag vect ~=Ofor Matlab-style vectorized 
% operations: faster, but requires extra temp store 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n=max(size(om)); % Number of modes considered 
w=zeros(2*n); % Initialization (saves time later) 

% 

if vect = 0 
% 

% "Standard" loop over matrix locations 
% 

% Compute each (2 x 2) Grammian block in turn. 


% and store in correct upper & lower locations in W 
% 

for i=l:n 
omi=om(i); 
iw=2*i-l; 
forj=i:n 
omj*om(j); 
jw»2*j-l; 
tl=ze*(omi+omj); 
t2=(omj*omj)-(omi*omi); 
t3»2*omi*omj; 
t4s=t3*tl; 

dij=2*t3*(tl *tl >+<t2*t2); 

wij*[t4 cobs*omj*t2; -cobs*omi*t2 14]; 

wij=K(b(i,:)*b<j,:)Vdij)*wij; 

w(iw:iw+l j w:j w+1 )=wij; 

w(jw:j w+14w:iw+l )=wij’; 

end; 

ifb(i,:)*b(i,:) , <eps % Avoid singularity 
w(iw:iw+l 4w:iw+l >=eps*eye(2); 
end; 
end; 
else 
% 

% Calculations in Matlab-style vectorized form: 

% first compute contribution from B 
% 

il=(l:2:2*n-l); 

i2=(2:2:2*n); 

% 

som=size(om); 

if som(l) >1 % om entered as a column vector 

tlsre^om^oDesCM}^^*^!)*©!^); 
t2=(ones(n,l )*(om. A 2)'-(om. A 2)*ones( 1 41 )); 
omj=ones(n,l)*om'; 
t3=2*om*om t ; 

else % om entered as a row vector 

t l=ze*(om , *one»( 1 ,n)+oDes(n, 1 )*om); 
t2=(ones(n,l)*(om A 2Hom. A 2)'*ones(l4i)); 
omj =ones(n, 1 )*om; 
t3=2*om’*om; 
end 
% 

d«2*t3*(tl A 2)+{l2 A 2); 

beta=b*b’; 

% 

% Avoid singularities caused by zero entries in b,c 
% 

betad=max(diag(beta), eps*ones(n,l)); 
beta=beU-diag(diag(beU)); % Zero all diagonals 
betasbeta+diag(betad); % Put back diag, or eps 
% 

w(il 41)*bcta.*t3.*tl 7d; 
w(il42)«cobs*beta.*omj.*t27d; 
w(i241)»-cobs*beta.*omj , .*t27d; 
w(i242>=beta.*t3.*tl 7d; 
end 
# 

####### 

# 

function [isort, nsub]=subsys(ro,romin); 
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% An M-function to determine those modes which are 
% deemed to be correlated, based on a given matrix 
% of modal correlation coefficients and a 
% specified threshold value 
% 

% Outputs: 

% 

% isort sorting required to produce subsystems 
% nsub: dimensions of the subsystems 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n=max(size(n>)); % Number of modes considered 
% 

% Perform various initializations 
% 

niim-0; % Total number of modes grouped 
nsubsQ; % No subsystem dimensions yet 
isortsfl:]; % Modal oidering: unchanged so far 
kvect*(n+ 1 )*ooes( 1 ji); % Makes unsorted modes last 
kflag-O; % 1 for subsystem 1, 2 for subsys 2, etc. 

% 

% Characterize each subsystem in turn 
% 

while ntim < n 
% 

% Initializations for this subsystem 
% 

nsubMfc % No modes yet found 

kflagskflag+1; % Increment flag 
itcst-iaort(num+l); % Mode to be tested 1st 
jtesriaort(nuin+l:n); % Test for coir’n toi 

% 

while iaempty(itest) 0 

% 

% itest contains a set of modes to be tested 
% 

inew«Q; % No modes for next pats yet 

% 

for i«itest 
for j»j test 

if ro(ij) >■ romin A kvectfi) > n 

% 

% Mode j is a new mode, correlated 

% toi: store in this subsystem 

% 

kvectQHcflag; 
nsubkansubk+1; 
if j i incwmfinew j]; end; 
end; 
end; 
end; 

% 

% Pick up new set of modes (if any) to test 
% 

itesUdnew; 

end 

% 

% This subsystem is finished: store its data 
% 


num*num+nrubk; % Total number of modes 
nsui^[nsub nsubk]; % Subsystem dimensions 
[kvsort isoit]*sort(kvect); % Sorted modes 
end 


function ro»occalc(w,vect); 

% 

% An M-function to compute the overall 
% correlation coefficients of a unifonnly- 
% damped flexible structure 
% 

% The flag vect -a 0 for Matlab-style vectorized 
% operations: Aster, but requires extra temp store 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n*max(size(w»/2; % Number of modes considered 
% 

if vect — 0 
% 

% "Standard* loop over matrix locations 
% 

ro“«yc(n); % Initialization (saves time later) 
forialm 
iw«2*i-l; 

wiiAnonn(w(iw:iw+l 4wriw+ l),'fro'); 
forj*i+l:n 
jw*2^j-l; 

wyfaaonnfwQwyw+ljwgw+lX’fro*); 
wijf=nonn(w(iw:iw+ 1 jw:j w+ 1 ),*fro f ); 
ro(iJ)*wijf/sqrt(wiif*wijf); 
ro(j>i)-io(ij ); 

end; 

end 

else 

% 

% FrObenius norm calculation in Matlab vector fonn 
% 

n2*2*n; 

e-eye(n2}+diag(oncs(n2-l,lX-l); 

e-c(:,l:2:ii2-lX 

% 

% Calc u late Ffobenius norms of each (2x2) block 
% 

ro«e , *(w.*w)*e; 

ro-sqrt(ro); 

% 

% Now normalize 
% 

rodiags-diag(aqrt(diag(io))); 

% 

ro-rodiags\roAodiags; 

end 

* 

# 

function 

[amodal,bmodal^modal]assmodal(on^ze,b,c); 

% 
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% An M-function to construct the modal state 
% space model corresponding to a uniformly- 
% damped flexible structure 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n=max(size(om)); % Number of modes considered 
% 

amodal=zen>s(2*n); 

adiag=[-2*ze -1;1 0]; % "Template" for diags of a 
% 

m=size(b); m=m(2); 
bmodal=zeros(2*n^n); 

% 

P=azc(c); p=p(l); 
cmodal=zeros(p,2*n); 

% 

for i=l:n 
i2=2*i-l; 

amodal(i2:i2+142:i2+l>=om(i)*adiag; 

hmodal(i2,:)=b(i,:); 

cmodal(:42>-c(:4); 

end 

# 

# 

function 

[am,bm,cm]=blkmult(ain,bm,cm,wctot4Wsort4isub); 

% 

% An M-function to apply a similarity 
% transformation stored as ordered blocks 
% to a given state space model 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

kmax=max(size(nsub)); % Number of subsystems 
% 

iwl=l; 
for k=l:kmax 

nsubk2»2*nsub(k); 
iw2=iw 1 +nsubk2- 1 ; 

iwvect=iwsort(iw 1 :iw2); % Required state order 

% 

% Retrieve k-th block of similarity transformation 
% 

wck=wctot( 1 :nsubk2,iwl :iw2); 

% 

% Premult row blocks of am, bm by inverse of wck 
% 

am( i wvect, : )=wck\am(i wvect, : ) ; 
bm(iwvect,: )=wck\bm(i wvect, : ) ; 

% 

% Postmultiply column blocks of am, cm by wck 
% 

am(:,iwvect>*am(:4wvect)*wck; 

cm(:,iwvect)-cm(: 4 wvect)*wck; 

% 

iwl*iw2+l; 

end 

# 


# 

function y=modstep(om*ze,b,c4u,t); 

% 

% An M-function to compute the step response 
% of a state space model in "symmetric" 

% modal form of a flexible structure 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

n=max(size(om)); % Number of modes 
tmax=max(size(t)); % Length of time vector 
p=size(c); p=p(l); % Number of outputs 
di=zeros(p,l); 

% 

y=zeros(tmax,p); 
adiag=[-2*ze -1; 1 0]; 

% 

for i=l:n 
% 

% Setup required submatrices 
% 

ai=om(i)*adiag; 
bi=[b<i 4 u); 0]; 
ci=[c(:4), 0*c(:4)]; 

% 

% Add step response of this mode to total 
% 

y=y+step(ai,bi,ci,di,l ,t); 
end 

# 

* 

function y=blkstep(a,b,c4u,t); 

% 

% An M-function to compute the step response 
% of a state space model of a flexible structure 
% with A block diagonal (modulo ordering!) 

% Trevor Williams, NASA JSC, August 19, 1992 
% 

tmax=max(size(t)); % Length of time vector 
p«size(c); p=p(l); % Number of outputs 

di=zen>s(p,l); 

% 

y=zeros(txnax*p); 

% 

% Determine OVERALL block structure of a, directly 
% 

[iasort4iasub]ssubsys(abs(a>4-ab8(a f ),eps); 

kmaxsmax(size(nasub)); 

% 

il-1; 

for k«l±max 

i2=il+nasub(k>-l; 

ia»iasort(il:i2); 

% 

% Add step response of this block to total 
% 

ysry+step(a(ia,ia),b(ia4uXc(:4a),di, 1 ,t); 

il«i2+l; 

end 
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